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1 Introduction 


Our purpose is to give mathematical studies for new 
algorithms applied to MEG. Concerning the direct 
problem, we follow the formulations of [1-5]; 
Maxwell's equation controls the total current density 
J(x) and the magnetic field B(x ) in such a way 
that V*5 = 0 and VxB = p 0 J. We assume that 
J = J p - o(x)VV , where J P (x) is the neuron current, 
and E = -VV stands for the electric field caused by 
J. The brain in consideration is denoted by a 
bounded domain £2 <z R 3 with boundary d£2. 


Therefore, J P (x) is a vector field on £1 with the null 
normal component on d£2. An important assumption 
is made on the conductivity g(x); g(x) = Gi (x e £1) 
and g(x) = g q (x e £2 c ), where Gi >0 are nonnegative 
constants. If £2 is a sphere and G 0 = 0, then we have 
a formula of [2] for J P (x) to determine B(x) directly. 
In the classical theory, it is assumed that the primary 
current J P (x) is a combination of dipoles such as 

Jp(x) = Yl k=l Q k 8{x-a k ). 


Unknown parameters {a k , Q k } are determined to 
adjust with discretely measured values of {n *B) on 
d£2 in use of the method of least square 
approximation ( LSA ). One aspect of the method of 
current element distribution is a trial to recover the 
total current density J P (x). First observation is that 
if J P (x) is a single dipole Q k S(x - a k ), then it 
produces the magnetic field 

M„ Qk x (x-a k ) 
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Let the number of channels, their positions, and the 
directions be M, xj e d£2, and Vj e S 2 , respectively, 
where 1 < j <M (usually, around 100 observing 
points are taken). Taking N sufficiently large 
(around 200), one introduces the mapping 
cp : ( Q, a) e R 6N ^ 
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and the functional 

&(Q,a) = ^\\(p(Q,a)-z\\ 2 RU , 

where 

Q = (QuQ 2 ,-QN)eR 3N , 

a = (a\, a 2 , ••• a N ) e R 3;V 
are parameters to be determined, and 

z = (z h z 2 , - z M ) e R M 

with zj g R standing for the observed datum of the v y 
component of the magnetic field B(x ) measured at 
the jth channel. Each Q k 8(x - a k ) is called the 
current element. First, N random current elements 
are distributed on £2 . One element is selected 
randomly and perturbed so that decreases. Then 
the next one is chosen also randomly to continue the 
process. If does not decrease, that element is left 
un-perturbed. If the directions of two adjacent 
elements happen to be opposite, then randomly 
selected one of them is moved away. This procedure 
is continued until a local minimum of $ is achieved. 

2 Second optimization 

Serious discussions have to be made on the relation 
between the number of assumed dipoles and that of 
channels; its balance leads LSA to be either under¬ 
determined or over-determined, of which difficulties 
are essentially different. 

Given cp : R” -» R w and z e R M , the discrete inverse 
problem is formulated as to find x e R” satisfying 
cp(x) = z, of which solution is called the strict 
solution. It is over-determined and under¬ 
determined according to n < m and n > m, 
respectively. In the over-determined case, the strict 
solution is hard to obtain; the discrete inverse 
problem is reformulated as the least square problem 
of finding xef satisfying fi(x) =y, where 

#(x) = ||<p(x) - z|| 2 and j = inf D. 



Its solution is called the least square solution , but 
practically one can examine its local minimality 
only. If fi(x) is satisfactorily small, one says that x 
has a high accuracy. Any local minimum with high 
accuracy is regarded as a least square solution. 
Sometimes it is not unique. One has to select a 
particular solution, regarding a priori information 
and their accuracies. The least square problem is 
said to be quasi-identifiable if any local minimum 
of with high accuracy is isolated. This allows one 
to select one of them as a desirable solution. A 
criterion for quasi-identiflability, given by [6] is the 
high-accuracy and the rank condition. The latter 
means that (p' (x 0 ) * cp' (x 0 ) ; R" —» R” is non¬ 
singular. Because the rank condition is generically 
holds by n < m, only high-accuracy is to be 
examined in that case. 

However, the current element distribution method is 
set to be under-determined. Generically, the set of 
quasi-solutions M = {x\ cp(x) = z } forms a manifold 
of n - m dimension. The iterative sequence {x/} 
approaches M. Then it freezes in the region of 

Ax /+1 ~ $(x) for Ax /+1 = x /+1 - x l , because in this 

range the accuracy is hard to improve. Making use 
of the tangent space T X M as the perturbation set, 

one can remove the freezed sequence without losing 
the high-accuracy. This process, called melting , 
is realized by the singular decomposition of the 
matrix (p ? (x/), because f>(x/) « 1 implies 

T Xl M ~ Ker<p' (x t ) . 

In the actual process of MEG, the above processes 
are realized in the following way. First, in 
approaching, mainly Q e K 3N is perturbed. Also, the 
maximum descendent method is better than the 
random perturbation here. On the other hand, in 
melting, mostly a e K 3N is moved. Furthermore, 
current elements are so taken to be clustered. We 
call this special process of melting the binding. One 
method for binding is as follows. First, for 8 > 0 
small, we take an 8 covering of 
S = e R 3 11 < k < a}. Fixing 8 > 0, this process 
is repeated in use of the random perturbation until 
the least number of the covering balls is obtained. 
We call it covering. Then, those balls are contracted 
as small as possible, remaining to be a covering E. 
We call this process the biting. After then, the 
melting is done under the constraint that any 
element does not move away from those contracted 
balls. Taking smaller 8 > 0, we repeat the process. 


For computations indicated by the next figures, a 
more rough method is adopted. There, only one ball 
is taken to contain E. Based on this covering, biting 
and melting are done similarly. Observing points are 
16, and current elements are 30. Then it was seen 
that current elements are clustering near the 
originally clustered 8 dipoles. 
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